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1. INTRODUCTION 

An intraguild predation (IGP) system refers to simultaneous predation and competition interactions 
where both IG predator and IG prey species utilize a common resource while also took benefit from preying 
upon one another [1]. The IGP phenomenon is common in nature and its important implications have been 
discussed to a large extent from both theoretical and practical perspectives [2]-[4]. IGP systems can exhibit 
complex and shifting dynamics between different contrasting behaviors which include among others the 
existence of multiple stable states (e.g. species extinction and coexistence) as well as random shifts among 
such different dynamics [6], [7]. As the number of species in an IGP system increases, the overall dynamical 
complexity and observed biodiversity also increase and eventually shape the functioning of an ecosystem [8]. 

To date, studies on the dynamics of IGP systems have attracted significant interests from researchers 
across various disciplines. Since its introduction in [9], various modeling and analyses of complex dynamics 
arising from IGP systems have been investigated in several studies, see for examples [[10J-[[14]. The work in 
demonstrates how IGP systems are found to be widspread in various ecosystems and significantly shape 
the trophic structure/biodiversity and the stability/persistence of an ecosystem. The work in reports that the 
resource preference of IGP predators play important roles in determining which predators occupy a particular 
ecosystem. More recently, an investigation of the impact that prey abundance and diversity have on the strength 
of predator-prey interaction in IGP systems are reported in [18]. Despite these developments, studies about IGP 
dynamics in the presence of alternative resources remain limited. 

This paper proposes a model of an IGP system with exclusive alternative resource and evaluates its 
qualitative dynamical properties using concepts from stability and bifurcation theories [19], [20]. Specifically, 
this paper identifies the parameter values that ensure the survival/coexistence of both IG predator prey in the 
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presence of alternative resources. The main finding in this paper is that if the IG predator is less efficient than 
the IG prey in converting the consumed shared resource into its growth utilizations, then the IG predator should 
satisfies a minimum predation rate parameter value ensure its survival and co-existence in the ecosystem. 


2. AN IGP SYSTEM MODEL 

Consider an asymmetric IGP system as shown in Figure [I] which consists of one predator and one 
prey which consume three resources. Let Nj and N2, respectively, be the IG prey and IG predator biomasses. 
Assume both species compete for a common resource No which is essential for their growth/reproduction. 
Each species N; and No also has an exclusive alternative resource N3 and N4, respectively, that is essential for 
survival. All resources satisfy logistic growth functions and their consumptions by the predator/prey species are 
modeled with type 2 response function [7], [21]. The consumption of the IG prey by IG predator is modeled as 
type 3 response function, while the death rates of both species are linearly proportional to the number/density 
of each species. The IGP system model considered in this paper is therefore defined as (ip. 


dNo/dt = No [ro (1 — (No/Ko)) — a10N1/(ho + No) — a20N2/(ho + No)| 

dN, /dt = N; [e10410No/(ho + No) — €13013N3/(h3 + N3) — c21NiNo/(q° + NÊ) — dı] 

dN2/dt = Nz [e20a20N0/(ho + No) — e24a24N4/ (h4 + Na) + c21 N? /(@? + NÊ) — do] (1) 
dN3/dt = N; [r3 (1 — (N3/K3)) — ai3Ni/(h3 + Na)| 

dN4/dt = N4[ra (1 — (Na/K4)) — az4N2/ (h4 + Na)| 


In (i). parameters r;, K;, and h; denote, respectively, the growth rate, carrying capacity, and half saturation 
constant of the ¿th resource (for i = 0,3, 4), while aj; and ej; denote the grazing rate and grazing efficiency of 
species 7 on resource j (for j = 1, 2), respectively. Furthermore, d; is the death rate of species j, and cx is the 
attack rate of species j on species k = 1,2 with k Æ j. The time variable is denoted with t. 


IG Predator 


IG Prey 
/ G TE 


Alternative Common Alternative 
Resource 1 Resource Resource 2 


Figure 1. Schematic of an asymmetric IGP system with an exclusive alternative resource 


Consider a non-dimensional parameterization [21], [22] of state variables dan parameters of (1) as (2). 


zo = No/Ko, rv, = M, tq = No, x3 = N3/Ks, v4 = N4/Ka, T = rot, 
ko = ho/Ko, ks = h3/Ks, k4 = h4a/Ka, 73 =13/T0, y4 = 14/10, c = C12/To, 
aio = @10/Koro, 20 = a20/Koro, ka = ha/Ka, Y3 = r3/To, ya = r4/T0, c= &12/T0, 
a= a10€10/To, €2 = a29€20/To, €3 = €13413/To, €4 = €24024/To.- 

(2) 


Define t; = dx,/dt. Then a parameterization of the state equation model (1) can be written as (). 


to = Xo [(1 — zo) — a1041/(ko + zo) — a20%2/(ko + X0)] 

i1 = q1 [e1x0/(ko + £0) + €3%3/(kg + £3) — cx £2/(q? + x3) — 5; | 

i2 = £2 |e2x£0/(ko + xo) — c4x4/(ka + £4) + 81 42/(q? + 27) — 59] (3) 
&3 = z3 [ys (1 — z3) — ai3@1/(k3 + z3)] 


t4 = x4 [y4 (1 — 24) — a24£2/ (k4 + x4)] 
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All parameters in are chosen to satisfy predation rate efficiency constraints. For instance, |e;| < 
|a;,;| must hold to account for conversion efficiency of species i in consuming resource j, or |ajo| < 1 to 
account for shared consumption of resource 0 by species 7. Assumptions below are also used in model 3). 

- Assumption |: the IG predator is less competitive than the IG prey in consuming the shared resource 


- Assumption 2: the IG predator is also less efficient than that of the IG prey in converting the consumed 
shared resource into their energy use or reproduction 
- Assumption 3: the IG predator and prey have similar consumption rate on their exclusive resource 
Assumptions 1 and 2 (as defined by setting €2 < €; ) suggest that IG prey may survives if its growth 
from consuming the shared resource is sufficient enough to overcome the predation pressure from IG predator. 
The IG predator, on the other hand, should maintains certain level of predation rate on IG prey to balance its 
loss from being less competitive in the consumption and less efficient in the conversion of the shared resource. 


3. METHOD 
3.1. Equilibrium points characterization 
To examine the stability of model (3), we first computed and identified its 13 possible equilibriums 
E; (i =1,...,13) which consist of 8 trivial and 5 nontrivial equilibriums as (4)-(6). 
a) Trivial equilibriums 


Eı = (0,0,0,0,0], £2 = [1,0,0,0,0], Æ = [0,0,0,1,0], £4 = [0,0,0,0, 1], 


E _ E _ (4) 
Es = (1,0, 0, 1, 0], Es = (1, 0, 0, 0, 1], Er = [0,0,0,1,1], Eg = [1,0,0, 1,1]. 

b) Nontrivial equilibriums 

= Eg = [z8 0,23, i3], where a = ko (ôi = (e3/kg + 1)) /[e1 = ôi + (e3/kg + 1)], 
x3 = (1 = x0) (ko + Xo) /Q20; and r3 = (1 = k4)/2 a ya = k4)? —4 (—k4 + 0249 /Ya). 

- Eio = [et 4,0,1., 1], where qao = [ko (ôi = (e4/ka + 1))]/lex = OT + (e4/k4 + 1)], 
xi = (1 = x0) (ko + xo) /Q20; and a = (1 = k3)/2 a ya k3)? 4( k3 } @24z1/73). 

= E, = ot at ar | where x}! = X (cf. for detailed formula), x}! = y3k3/a13, x}! = 
—ky + Q24X2/Ya(1 = ta) and xil => (6 = €xo/(ko + zo)) /cy3k3/013 (¢? + (y3k3/013)) ‘ 


- Ei = [x§?, 21°, x57, 73", 0], where «3? = y4k4/a24 while the other elements satisfy (5). 


2 2 
e2xg /(£0° + ko) +c (ai) /[a? + (21°) ] -— & = 0, 
xy” — (ko + 297) (1 — 297) /a19 — Y4k4a20/a24010 = 0, (5) 
(kg + ©3”)(1 — 23°) — aisr] /73 = 0. 


- Eis = (xp°, x73, 75°, 133, 74°], where each equilibrium’s element satisfies (6). 


(1 — 23) — aorist? / (ko + x°) — a20%025° 247 / (ko + rozi?) = 0, 
exo /(ko + x0) — €323°/(k3 + 13°) — exp?a3?/(¢° + (x7°)?) — ô = 0, 

egg’ /(ko + xq) — eaz3 /(ka + 24°) — eler)” /(@ + (@7°)”) — d2 = 0, (6) 
y3(1 — x3?) — ai3ay°/(kg + 23°) = 0, 


ya(1 — rI’) = aoar} /(ka + zi’) =0. 


3.2. Local stability analysis 
Based on the equilibriums data obtained in section 3.1, the system’s local stability property at each 
equilibrium can be analyzed using linearization method. The main objective in this analysis is to examine 
which equilibrium points will ensure stable species co-existence of both IG predator and IG prey species. 
Consider nonlinear system model in (3) which satisfies the following general form. 


tolt) = fo(zo, peas , In) 
(7) 


En(t) = Fr{Zo; tee , Tn), 
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The linearization of (7) around each of its equilibrium point Ey (£ = 1,--- , 13) is defined as (8). 


Z(t) Sho vee oe £o(t) Ai, Ain Z(t) 

: = 3 . 3 — | 3 h; l : = Az(t) (8) 
> Ofn afn = z 
Ey(t) Sta a fa ueis Tn(t) Ani ce Ann (=E) Talt) 


The elements of the linearized system matrix A in are detailed in [23]. The local stability property of 
the IGP system at each equilibrium point may then be evaluated using (8), and the results are summarized in 
Table [I] It can be seen in this table that most of the equilibrium points of model (3) are unstable. 


Table 1. Local stability properties of equilibrium points of model 


Equilibrium Eigenvalues Stability 
Eı = (0, 0,0, 0, 0]7 —0.0100, —0.0100, 1.0000, 1.0000, 1.0000) unstable 
E> = [1,0,0, 0, 0] 0, —1.0000, 0.0138, 0.0138, 1.0000, 1.0000) unstable 
E3 = [0,0,0,1,0]" 0, 1.0000, —1.0000, —0.0094, —0.0100, 1.0000) unstable 
E, = [0,0,0, 0, 1)” 1.0000, —0.0100, — 1.0000, —0.0094, 1.0000) unstable 
Es = [1,0,0, 1, 0]7 —1.0000, —1.0000, 0.0144, 0.0138, 1.0000) unstable 
Es = [1,0,0, 0, 1] —1.0000, 0.0138, —1.0000, 0.0144, 1.0000) unstable 
E7 = [0,0,0, 1, 1)" 1.0000, — 1.0000, —0.0094, — 1.0000, — 0.0094) stable 
Eg = [1,0,0,1, 1)" —1.0000, —1.0000, 0.0144, —1.0000, 0.0144) unstable 
Eg 0.0164, —0.3546, —0.6856, — 1.0000, 0.0000) marginally stable 
Eo —0.6614, 0.0172, —492.0332, —1.0000, 0.0100) unstable 
Fur —6.8750, 0.7178, —0.0087, —0.0047, 0) unstable 


4. RESULTS AND DISCUSSION 

Since the local stability analysis for a particular parameter set suggests most equilibrium points to be 
unstable, bifurcation analysis is further conducted to search for possible existence of multiple stable equilib- 
rium points for different combinations of parameter values. The existence of such multiple (nonzero) stable 
equilibriums essentially indicates the existence several states where all species can survive and co-exist. 

The computation of stable co-existence equilibriums of (3) is challenging as it requires the solution 
of polynomial function equations of order 6 with 18 parameters. To compute possible co-existence equilib- 
riums, bifurcation analysis which starts at a nominal, stable species co-existence equilibrium was conducted. 
The following initial parameter values are chosen in the bifurcation analysis: kọ = 20, k3 = k4 = 15, 
a0 = 0.6, Q20 = 0.4, 13 = 0.1, a24 >= 0.15, c= 0.01, qd = 2, Oy = 69 = 0.01, yg = jf4 = 1, 
€ = 0.5, €2 = eg = 0.1, and e3 = 0.01. Figure [2] shows the resulting trajectories and phase portraits, es- 
pecially in Figures [2ha) and Plo), Figure P| show a stable co-existence equilibrium at [x£o, £1, £2, £3, £4]? = 
[0.8545, 2.1211, 4.4054, 0.9867, 0.9724]* is reached. This stable point is therefore used as the starting point 
for the one-parameter bifurcation analysis of model (). The bifurcation analysis reported in this paper were 
obtained using MATCONT toolbox in MATLAB [25]. 
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55 


Figure 2. The IG predator’s (a) time response and (b) phase portrait, the IG prey dynamics for the nominal 
model parameter set with initial condition [x9, x9}, z8, z8, x9]? = [1,5,5,0.1,0.1]”) 
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Note that a bifurcation point is essentially a point where the number and qualitative behaviors of an 
equilibrium change as the parameters are varied. A Hopf (H) bifurcation is a point where an equilibrium 
changes from stable to unstable, or vice versa. A Limit Point (LP) bifurcation is a threshold point in which an 
equilibrium point jumps between two stable states with contrasting qualitative behavior. Finally, a Branching 
Point (BP) is the point from which periodic doubling of equilibriums (such as in logistic map) starts to occur. 


4.1. Bifurcation analysis to identify species coexistence equilibriums 

The parameters that are examined in this analysis case are c; (ii = 1,...,4), di2 (i2 = 1,2), 
and c. Other parameters are kept constant. The results of one-parameter bifucation analysis are plotted 
in Figure The indicated INIT point in Figure Blis the initial equilibrium point with response shown in 
Figure|2} Each red-star mark is a bifurcation point whose corresponding symbol 0% indicates the i-th bifur- 
cation point of type 0 € {H, BP, LP} under variation of one bifurcation parameter j € {€;,6;,c}. Each 
bifurcation point is obtained by computing the forward (increasing) and backward (decreasing) variations of 
one bifurcation parameter while keeping the other fixed. For example, Figure[3]shows that forward and back- 
ward variations of parameter 6;2 produce three bifurcation points which include two Hopf (H}, H3) and one 
LimitPoint (LP) bifurcation points. Moreover, there exist multiple species co-existence attractors among 
which the IGP system may jump from one to another when one bifurcation parameter is varied. 

Table |2| summarizes all bifurcation points and their equilibrium points x* which occur in the one- 
parameter bifurcation analysis. Figure [4]illustrates the response gla) and phase portrait (4[b)) of model (3) at 
the bifurcation parameter ¢; = 6.0507 which differ significantly from the nominal response shown in Figure[2| 
The biological interpretation of the results in Figure[3]and Table[2]can be examined based on the corresponding 
bifurcation parameter. For instance, consider the impacts of variation on IGP prey’s consumption rate £; from 
its nominal value of £, = 0.5 at the indicated INIT point in Figure[3| It can be seen that forward computation of 
€1 drives the IGP system towards bifurcation points H, os and H, 2 which correspond to €; values of €; = 1.199 
and ¢, = 6.0507, respectively. However, from biological system modeling perspectives, these values are not 
reasonable as the consumption rate parameters should always satisfy a constraint of the form |e1| < 1. Thus, 
one may concludes that consumption rate parameter £; that is higher than the nominal value will not change 
the stability of IGP system (3p. Similar interpretation may also be inferred for other parameters. 


35 


Figure 3. One-parameter bifurcation analysis for species co-existence 


Table 2. Bifurcation points and the corresponding equilibrium points 


Parameter xro xi r3 x3 zi Bifurcation type 
€1 = 2.049 0.947785 1.822995 000021 0.988598 0.000011 Branching 
€1 = 1.199 0.690866 2.363633 12.4451 70.985214 0.882463 Hopf 
e1 = 6.0507 0.240142 4.199785 32.149418 0.973708 0.692697 Hopf 
€2 = 0.079 0.845958 2.395845 4.434099 0.985012 0.958322 Hopf 
e3 = 0.267 0.745452 2.234137 9.850565 0.986024 0.907112 Hopf 
61 = 0.0110003 0.850746 2.459088 4.091537 0.984616 0.961549 Hopf 
61 = 0.13525 0.763841 5.804019 3.552916 0.963642 0.966622 Hopf 
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Figure 4. Contact time effect on metal ions adsorption in (a) pseudo-second order kinetic model and 
(b) (M?+] = 0.14 mmol L-1, pH=6 at RT) 


4.2. Bifurcation analysis in the presence of consumer’s alternative resource 

One main question on the IGP model is whether the IG predator can survive by only consum- 
ing the shared alternative resource and without preying upon the IG prey. We examine this question for the 
case when the IG predator is less efficient than the IG prey in converting the consumed shared resource into en- 
ergy/reproduction uses. Formally, this case can be defined as a constraint on parameters of the form: |e2| < |e1l. 

The backward computation on predation rate c was conducted and the bifurcation curve of the IG 
predator density with respect to parameter c is plotted in Figure|5| This plot shows that, depending on the 
values of c, two types of equilibriums for the IG predator may occur, namely a stable equlibrium (horizontal 
curve) and an unstable equilibrium (vertical curve). On one hand, the unstable equilibrium curve shown in 
Figure [5] tends to result in IG predator extinction at a BP bifurcation point for the value of c = 0.007238. 
This suggests that c = 0.007238 is a lower bound value of predation rate c which should be maintained by IG 
predator to avoid extinction. On the other hand, the LP bifurcation which occurs at c = 0.006269 corresponds 
to the threshold value c which will guarantee the IG predator survival for some time before eventually extincts. 
Based on the results shown in Figure|5| one may thus defines the following three partitioned values of parameter 
c which also correspond to different equilibriums/behaviors of the IGP system dynamics: 
- if0 < c < 0.006269 (LP), then the IG predator will extinct without any chance to grow/reproduce. 


- if (LP) 0.006269 < c < 0.007238 (BP), then the IG predator will exist and survive for some period of 
time before eventually extincts. 


- if c > 0.007238 (BP), then the IG predator will survive and co-exist with other species. Its dynamics will 
exhibit limit cycles below c = 0.0085 and become asymptotically stable above c = 0.0085. 

From biological perspective, the partitioned values of parameter c suggest that, even when alternative 
resources exist, if the IG predator is less efficient than the IGP prey in converting the consumed shared resource 
into energy/reproduction uses, there still exists a lower bound value of c that the IG predator should maintains 
to ensure its co-existence survival. Moreover, even when the IGP predator consumption rate on alternative 
resource is c = 1 (100% efficiency), it cannot survive without predating on the IGP prey. One may thus 
concludes that, to ensure species co-existence in IGP system, the IG predator predation rate on IG prey is a 
more dominant factor than its consumption rate on the shared alternative resource. 

In contrast to the IG predator species, the IG prey consumption rate £3 on its alternative resource does 
not play significant roles to ensure its co-existence in the system. As shown in Figure [6] the IG prey density 
which starts at the INIT point will increase towards the LP bifurcation point as £3 is increased. However, since 
the constraint 0 < £3 < 1 must holds, the LP bifurcation point (€3 = 1.109 ) will never be realized. Thus, the 
IG predator is guaranteed to survive even without consuming the alternative resource. Specifically, the Hopf 
bifurcation point H4 where IG prey extincts will never be reached since the system dynamics must follow the 
bifurcation curve by first visiting points LP, H2, H3, BP which, in particular, are also not reachable. 
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Figure 5. Bifurcation curve for variation on parameter c 
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Figure 6. Bifurcation curve for variation on parameter £3 


5. CONCLUSION 

This paper has developed and analyzed the qualitative properties of an IGP system model in the pres- 
ence of exclusive alternative resource. Local stability and bifurcation analyses were performed on the model 
to examine the values of parameter set which guarantee the existence of stable equilibriums where both IG 
predator and IG prey species can survive and co-exist. One main finding in this paper suggests that, in the 
presence of alternative exclusive resource, if the IG predator is less efficient than the IGP prey in converting 
the consumed shared resource into energy/reproduction utilizations, then there is a lower bound on the value 
of predation rate parameter which the IG predator should maintains to ensure the co-existing survival of both 
the IG predator and IG prey in the ecosystem. The results presented in this paper also suggest that in order to 
ensure species co-existence in the ecosystem, the IG predator’s predation rate on IG prey is a more dominant 
parameter/factor than its consumption rate parameter on the alternative resource. 
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